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Abstract 

We  present  implementation  details  and  analyze  convergence  of  a  two-grid  solver 
forming  the  core  of  a  fully  automatic  hp-adaptive  strategy  for  electromagnetic  prob¬ 
lems  [9,  26].  The  solver  delivers  a  solution  for  a  fine  grid  obtained  from  an  arbitrary 
coarse  hp  grid  by  a  global  hp- refinement.  The  classical  V-cycle  algorithm  combines  an 
overlapping  Block-Jacobi  smoother  with  optimal  relaxation,  and  a  direct  solve  on  the 
coarse  grid.  A  theoretical  analysis  of  the  two  grid  solver  is  illustrated  with  numerical 
experiments.  Several  electromagnetic  applications  show  the  efficiency  of  combining 
the  fully  automatic  /ip-adaptive  strategy  with  the  two  grid  solver.  This  paper  is  a 
continuation  of  [25]. 
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1  Introduction 


The  paper  is  concerned  with  a  construction  and  study  of  an  iterative  solver  for  linear  systems 
resulting  from  ftp-adaptive  Finite  Element  (FE)  discretizations  of  Maxwell’s  equations.  Here 
ft  stands  for  element  size,  and  p  denotes  element  order  of  approximation,  both  varying  locally 
throughout  the  mesh. 

The  algorithm  presented  in  [13,  9]  produces  a  sequence  of  optimally  /ip-refined  meshes 
that  delivers  exponential  convergence  rates  in  terms  of  the  FE  error  measured  in  energy  norm 
vs  the  discrete  problem  size  (number  of  degrees-of-freedom  (d.o.f.))  or  the  CPU  time.  A 
given  (coarse)  hp  mesh  is  first  refined  globally  in  both  h  and  p  to  yield  a  fine  mesh ,  i.e.  each 
element  is  broken  into  four  element-sons  (eight  in  3D),  and  the  discretization  order  p  is  raised 
uniformly  by  one.  We  solve  then  the  problem  of  interest  on  the  fine  mesh.  The  next  optimal 
coarse  mesh  is  determined  by  minimizing  the  projection  based  interpolation  error  of  the  fine 
mesh  solution  with  respect  to  the  optimally  refined  coarse  mesh.  The  algorithm  is  very 
general,  and  it  applies  to  H1-,  H(curl )-,  and  H( div) -conforming  discretizations  [12,  10].  In 
particular,  it  is  suitable  for  electromagnetic  problems.  Moreover,  since  the  mesh  optimization 
process  is  based  on  minimizing  the  interpolation  error  rather  than  the  residual,  the  algorithm 
is  problem  independent,  and  it  can  be  applied  to  nonlinear  and  eigenvalue  problems  as  well. 

Critical  to  the  success  of  the  proposed  adaptive  strategy  is  the  solution  of  the  fine  grid 
problem.  Typically,  in  3D,  the  global  ftp-refinement  increases  the  problem  size  at  least  by  one 
order  of  magnitude,  making  the  use  of  an  iterative  solver  inevitable.  With  a  multigrid  solver 
in  mind,  we  choose  to  implement  first  a  two  grid  solver  based  on  the  interaction  between  the 
coarse  and  fine  hp  meshes.  The  choice  is  quite  natural.  The  coarse  meshes  are  minimum  in 
size.  Also,  for  wave  propagation  problems  in  the  frequency  domain,  the  size  of  the  coarsest 
mesh  in  the  multigrid  algorithm  is  limited  by  the  condition  that  the  mesh  has  to  resolve  all 
eigenvalues  below  the  frequency  of  interest.  Consequently,  the  sequence  of  multigrid  meshes 
may  be  limited  to  just  a  few  meshes  only. 

The  fine  mesh  is  obtained  from  the  coarse  mesh  by  the  global  ftp-refinement.  This  guar¬ 
antees  that  the  corresponding  FE  spaces  are  nested,  and  allows  for  the  standard  construction 
of  the  prolongation  and  restriction  operators.  Notice  that  the  sequence  of  optimal  coarse  ftp 
meshes  produced  by  the  self-adaptive  algorithm  discussed  above  is  not  nested.  The  coarse 
meshes  are  highly  non-uniform,  both  in  element  size  ft.  and  order  of  approximation  p,  and 
they  frequently  include  anisotropically  refined  elements  (construction  of  multigrid  algorithms 
for  such  anisotropically  refined  meshes  is  sometimes  difficult),  but  the  global  refinement  fa¬ 
cilitates  greatly  the  logic  of  implementation. 

Customarily,  any  work  on  iterative  methods  starts  with  self-adjoint  and  positive  definite 
problems,  and  this  was  the  subject  of  the  work  presented  in  [25].  We  included  2D  and  3D 
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examples  of  problems  with  highly  non  homogeneous  and  anisotropic  material  data,  as  well 
as  problems  presenting  corners  and  edge  singularities. 

In  this  paper,  we  are  concerned  with  a  construction  of  a  similar  but  yet  different  two 
grid  solver  algorithm  suitable  for  general  electromagnetic  (EM)  problems.  We  also  discuss 
advantages  and  limitations  of  the  fop-adaptive  strategy  combined  with  the  two  grid  solver 
when  applied  to  real  life  EM  problems. 

The  structure  of  our  presentation  is  as  follows.  We  begin  with  a  formulation  of  the 
two  grid  solver  algorithm,  and  a  study  of  its  convergence.  In  Section  3,  we  present  some 
implementation  details,  while  Section  4  is  devoted  to  numerical  experimentation.  A  number 
of  EM  applications  is  presented  in  Section  5.  Conclusions  are  drawn  in  Section  6. 

Notice  that  the  two  grid  solver  is  not  intended  to  be  used  only  as  a  solver  itself,  but  also 
as  a  crucial  part  of  the  fop-adaptive  strategy.  Among  several  implementation  and  theoretical 
issues  that  we  address  in  this  paper,  one  is  especially  important  for  us;  is  it  possible  to  guide 
the  optimal  hp-refinements  for  EM  problems  with  a  partially  converged  fine  grid  solution 
only,  and  to  what  extent? 


Here,  u  is  an  angular  frequency,  e,  /i,  a  denote  dielectric  permittivity,  magnetic  permeability, 
and  conductivity  of  the  medium,  J"np  is  a  prescribed,  impressed  (source)  current,  Jgnp  is  a 
prescribed,  impressed  surface  current  tangent  to  boundary  r2,  n-  Jlfnp  =  0,  with  n  denoting 
the  normal  outward  unit  vector  to  T.  Finally,  j  is  the  imaginary  unit. 

For  the  sake  of  simplicity,  we  shall  restrict  ourselves  to  simply  connected  domains  72  only, 
avoiding  the  technical  issues  associated  to  cohomology  spaces,  see  e.g.  [7]. 


Standard  variational  formulation.  The  standard  variational  formulation  is  obtained 
by  multiplying  (2.1)  by  a  vector  test  function  F,  integrating  over  domain  72,  integrating  by 
parts,  and  using  the  Neumann  boundary  condition. 

(  Find  E  G  77D(curl;72)  such  that 

/  —(V  x  E)  ■  (V  x  F)dx  —  [  (c u2e  —  jua)E  ■  Fdx  =  (2  41 

Jn  p  Jn  \  ■  J 

— ju  l  Jimp  ■  Fdx  +  ju  [  JTP  ■  FdS  for  all  F  G  HD (curl;  O)  . 

Jn  J  r2 

In  the  above  /7z)(curl;  f2)  is  the  Hilbert  space  of  admissible  solutions, 

Hd( curl;  O)  :=  {E  G  L2(Q)  :  V  x  E  G  L2(n),  n  x  E  =  0  on  Tj  ,  (2.5) 

with  inner  product  defined  by: 

(u,  v)hd( curi;a)  :=  («,  v)L2(n)  +  (' V  x  u.  V  x  v)L2(Q)  .  (2.6) 

The  original  and  variational  formulations  are  equivalent  to  each  other. 


Stabilized  variational  formulation.  Introducing  a  space  of  Lagrange  multipliers  (scalar 
potentials): 

H^Ll)  :=  {q  G  H\n)  :  q  =  0  on  T, },  (2.7) 


we  employ  a  special  test  function  F  =  Vg,g  G  77^(11),  to  learn  that  solution  E  to  the 
variational  formulation  must  automatically  satisfy  the  weak  form  of  the  continuity  equation, 

-  /  (oA  -  jua)E  ■ Vqdx  =  -ju  t  Jimp  ■ Vqdx  +  ju  \  jfp  ■  VqdS  .  (2.8) 

Jn  Jn  Jr2 

We  also  recall  the  Helmholtz  decomposition: 

E  =  V0  +  E0,  where  (j)  G  77^(0)  and  (E0,  ~Vq)  =  0  Wq  G  77^(71)  .  (2.9) 
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It  is  well  known  that  the  standard  variational  formulation  is  not  uniformly  stable  with 
respect  to  the  wave  number  k 2  =  p(u2e  —  juxr).  As  k  — >  0,  we  loose  the  control  over 
gradients.  This  corresponds  to  the  fact  that,  in  the  limiting  case  k  —  0,  the  problem  is 
ill-posed  as  the  gradient  component  remains  undetermined.  A  remedy  to  this  problem  is  to 
enforce  the  continuity  equation  explicitly  at  the  expense  of  introducing  a  Lagrange  multiplier 
p  G  Hln(Q).  The  so  called  stabilized  variational  formulation  looks  as  follows. 

Find  E  G  iLi)(curl;  £l),p  G  HjjfD.)  such  that 

/  —(V  x£)(V  x  F)dx  —  [  (ca2e  —  jua)E  ■  Fdx 

Jn  p  Jn 

-  f  (cA  -  jua)Vp  ■  Fdx  =  -ju  f  Fmp  ■  Fdx  +  ju  f  jfp  ■  FdS 

<  J o  Jn  J r2  (2.10) 

VF  G  Hd( curl;Q) 

-  /  (cue  -  ja)E  ■ Vqdx  =  -j  f  J,mp  ■ Vqdx  +  j  f  jfp  ■  VqdS 

J  J  jt2 

Mq  G  HlD{9)  . 


By  repeating  the  trick  with  the  substitution  F  =  'Vq  in  the  first  equation,  we  learn  that 
the  Lagrange  multiplier  p  identically  vanishes ,  and  for  that  reason,  it  is  frequently  called  the 
hidden  variable.  In  comparison  with  the  original  formulation,  the  stability  constant  for  the 
regularized  formulation  converges  to  zero  slower  as  k  — >  0.  In  the  case  when  cr  =  0,  and 
the  right-hand  side  of  the  second  equation  vanishes,  we  can  divide  the  second  equation  by 
another  u  to  obtain 


/  eE  ■  Vg  dx 

Jn 


0. 


(2.11) 


In  this  case,  the  inf-sup  stability  constant  converges  to  one  as  u>  — >  0.  The  regularized  for¬ 
mulation  works  because  gradients  of  the  scalar- valued  potentials  from  Hxn{ii)  form  precisely 
the  null  space  of  the  curl-curl  operator. 

The  point  about  the  stabilized  (mixed)  formulation  is  that,  whether  we  use  it  or  not  in 
the  actual  computations  (the  improved  stability  is  one  good  reason  to  do  it),  the  original 
variational  problem  is  equivalent  to  the  mixed  problem. 


2.2  Formulation  of  the  two  grid  solver  for  EM 

Solving  a  linear  system  of  equations  (using  a  multigrid  scheme)  arising  from  Maxwell’s 
equations  is  challenging  mainly  for  two  reasons:  the  linear  system  is  (in  general)  indefinite, 
and  the  null  space  of  the  differential  operator  curl  is  large. 
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The  problem  of  indefiniteness  of  the  linear  system  can  be  overcome  by  requesting  the 
coarse  grid  to  be  fine  enough  (see,  for  example,  [6,  16]).  This  assumption  is  needed  both  to 
define  a  block-Jacobi  smoother,  as  well  as  to  prove  convergence  of  the  overall  two  grid  solver 
algorithm. 

In  order  to  control  the  solution  over  the  null  space  of  the  curl,  we  may  utilize  Helmholtz 
decomposition  (2.12),  and  treat  both  terms  separately. 

curl;  O)  =  (A'er(curl))  ©  (K er(curl))_L  (2-12) 


Two  corresponding  decompositions  have  been  constructed  for  lower  order  FE  spaces. 
More  precisely,  let  T  be  a  grid,  M  the  associated  lowest  order  Nedelec  subspaces  of  F/£>(curl;  0) 
of  the  first  kind  [24],  and  W  the  corresponding  first  order  piecewise  polynomial  subspace  of 
H\{Lt).  Let,  Vi  (resp.  eA  denote  the  non-Dirichlet  vertexes  (resp.  edges)  of  the  grid  T. 


Then,  we  define: 

Of  =  int(lJ{L  e  T  :  vt  e  dL }) 

Of  =  int(|J{Z  e  T  :  et  e  dL})  . 

And: 

Mi  =  {u  G  M  :  supp(w)  C  Of}  ;  Mf  = 
Wl’  =  {ueW  :  supp (u)  C  Of}  ;  W?  = 


(2.13) 

(2.14) 

M  : supp(u)  C  Of} 

(2.15) 

W  :  supp(w)  C  Of}  =  0  . 

(2.16) 

We  have  the  decomposition: 

W  =  Y.W’-  (2.17) 

V 

Arnold  et.  al.  [2]  proposed  the  following  decomposition  of  M: 

M  =  J2Mi  ,  (2-18) 

V 

which  we  shall  call  the  AFW  decomposition.  Another  well  known  decomposition  of  M  is 
Hiptmair’s  decomposition  [19]: 

M  =  £  M?  +  £  vwr .  (2.19) 

e  v 
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Each  decomposition,  together  with  the  already  prescribed  coarse  grid,  determines  a  two 
grid  solver  in  terms  of  a  multigrid  framework,  as  presented,  for  example,  in  [5].  More 
precisely,  the  bilinear  form  defined  over  each  subspace  can  be  inverted,  generating  a  block 
Jacobi  (or  Gauss-Seidel)  smoother  for  the  fine  grid  that,  together  with  the  coarse  grid,  define 
a  two  grid  solver  algorithm. 

A  formal  generalization  of  these  decompositions  for  /?,p-edge  elements  is  straightfor¬ 
ward.  Notice  that  Hiptmair’s  decomposition  (with  lowest  order  elements)  utilizes  only 
1-dimensional  subspaces  (and  therefore,  point  smoothers),  while  the  AFW  decomposition 
utilizes  4-dimensional  subspaces.  For  higher  order  elements,  size  of  patches  will  become 
considerably  larger  as  p  increases  and,  as  a  consequence,  amount  of  memory  (and  number 
of  operations)  required  by  the  corresponding  block  Jacobi  smoothers  become  prohibitive. 
Thus,  a  suitable  two  grid  solver  algorithm  for  hp-e dge  FE  may  come  from  combining  the 
ideas  presented  in  [25]  with  Hiptmair’s  approach  to  control  gradients. 

In  the  remainder  of  this  section,  we  present  a  two  grid  solver  algorithm  that  combines 
a  generalization  of  Hiptmair’s  approach  to  hp-e  dge  FE  with  the  block  Jacobi  smoother  pre¬ 
sented  in  [25].  M  and  W  will  denote  the  hp- FE  subspaces  of  //^(curl)  and  H^,  respectively. 

We  will  illustrate  via  numerical  experiments  the  importance  of  an  adequate  control  of 
the  kernel  of  the  curl  operator  formed  by  gradients  of  potentials.  Indeed,  a  two  grid  solver 
may  not  converge  if  the  gradients  are  not  resolved  correctly. 

Overlapping  block  Jacobi  smoothers.  At  this  point,  we  define  two  overlapping  block 
Jacobi  smoothers: 

•  one  used  as  a  preconditioner  for  the  electric  field,  given  by 

N 

sE  =  E  uopiT, 

l=i 

•  and  another  used  as  a  preconditioner  for  the  gradients,  given  by 

N 

sv  — 

i=i 

Here  Di  denotes  the  diagonal  sub-block  of  global  stiffness  matrix  A  corresponding 
to  d.o.f.  of  a  particular  (modified)  element  that  span  an  /?,p-edge  FE  subspace  Mi  C 
Hd(  curl;fl).  DVtl  denotes  the  diagonal  sub-block  of  global  mass  matrix  for  the  gradi¬ 
ents  (Laplace  equation)  Ay  corresponding  to  d.o.f.  of  a  particular  (modified)  element  that 
span  an  hp  FE  subspace  HJ  C  iJ^(fl).  are  the  matrices  associated  with  the  natural 

embeddings  from  Mi  into  M,  and  VHJ  into  M,  respectively. 
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At  this  point,  we  would  like  to  simplify  our  notation  and  drop  the  boldface  symbols  for 
the  matrices  and  vectors  in  the  coefficient  space. 


An  approximation  to  the  optimal  relaxation  parameter.  An  optimal  relaxation  pa¬ 
rameter  was  selected  in  [25]  to  minimize  the  error  in  the  energy  norm,  which  turned  out  to  be 
a  computable  number  for  self-adjoint  positive  definite  (SPD)  problems.  For  electrodynamic 
problems,  computation  of  the  optimal  relaxation  parameter  involves  solution  of  the  system 
of  linear  equations  that  we  are  trying  to  solve.  Thus,  only  an  approximation  to  it  may  be 
available. 

Since  S  ~  A-1,  we  define  our  approximation  to  the  optimal  relaxation  parameter  as  the 
argument  that  minimizes  Sr  in  the  energy  norm.  Thus,  at  step  n,  cd")  is  given  by, 


a 


(n) 


arg  min  ||JSV,(ri+1) 


(SAn\SASr^)B 
{SASr^),SASr(n))B  ’ 


(2.20) 


where  B  is  the  mass  matrix  associated  with  the  energy  norm  (( Bu,v )  =  (u,v)hd( curl)),  and 
S  is  either  Se  or  S y. 


2.3  The  two  grid  algorithm 

We  define  our  two-grid  solver  (based  on  a  modification  of  Hiptmair’s  decomposition)  as  the 
iteration  along  the  following  steps.  Given  current  solution  x,  and  residual  r,  we  perform, 

1.  coarse  grid  correction,  i.e., 

•  restrict  the  residual  to  the  coarse  grid  dual  space, 

r o  =  QTr  ;  (2.21) 

•  solve  the  coarse  grid  problem  for  coarse  grid  correction  Aa;o, 

AoAxo  —  r0  ;  (2.22) 

•  prolong  the  coarse  grid  correction  to  the  fine  grid  space,  and  compute  the  corre¬ 
sponding  correction  for  the  residual, 

Ax  =  QAxo,  A r  =  AAx  ,  (2.23) 

•  update  the  fine  grid  solution  and  residual, 


x  =  x  +  Ax  , 
r  —  r  —  A  r  ; 


(2.24) 


2.  block- Jacobi  smoother  on  the  fine  grid,  i.e., 


compute  the  smoothing  correction  and  the  corresponding  correction  for  the  resid¬ 
ual, 


Ax  =  Syr,  A  r  =  AAx] 

compute  an  approximation  of  the  optimal  relaxation  parameter  a, 

_  (Set,SeASet)b 
(■ SeASet  ,  SeASet  )b  ’ 
update  the  solution  and  residual, 
x  =  x  +  agAr  , 
r  =  r  —  «gAr  ; 


(2.25) 


(2.26) 


(2.27) 


3.  block-Jacobi  smoother  to  control  gradients,  i.e., 

•  compute  the  smoothing  correction  and  the  corresponding  correction  for  the  gra¬ 
dients  of  the  residual, 


Ax  =  S^r,  A  r  =  AAx] 

compute  an  approximation  of  the  optimal  relaxation  parameter  a, 

=  (Syr,  SyASyr)B 
(SyASyr,  Sy ASyr ) b  ’ 
update  the  solution  and  residual, 
x  =  x  +  ay  Ax  , 
r  —  r  —  ay  A  r  . 


(2.28) 


(2.29) 


(2.30) 


2.4  Convergence  Theory 

A  proof  of  convergence  for  our  two  grid  solver  algorithm  can  be  inferred  from  the  convergence 
theory  for  multigrid  algorithms  presented  in  [16],  which  refers  to  [15,  19],  and  [2]  among 
others,  for  detailed  proofs.  In  here,  we  outline  the  main  ingredients  of  the  convergence 
proof,  which  can  be  divided  into  three  parts. 

•  First,  we  introduce  some  notation  and  a  discrete  Poincare-Friedrichs  type  inequality, 
necessary  to  define  our  block  Jacobi  smoothers. 
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•  Next,  we  define  an  auxiliary  problem,  which  differs  from  our  original  problem  in  the 
value  of  wave  number  (squared)  k2.  By  setting  k2  =  —1,  we  obtain  a  SPD  auxiliary 
problem  with  convergence  properties  in  terms  of  the  two  grid  solver  equivalent  (up  to 
a  constant  times  element  size  h)  to  our  original  problem,  under  the  assumption  that 
the  coarse  grid  is  fine  enough. 

•  Finally,  we  prove  convergence  of  our  two  grid  algorithm  for  the  auxiliary  SPD  problem 
with  a  contraction  constant  independent  of  h,  and  possibly  depending  upon  polynomial 
order  of  approximation  p.  Thus,  convergence  of  the  two  grid  solver  for  the  original 
problem  is  guaranteed. 

We  assume  that  k2  is  real  and  such  that  our  boundary  value  problem  given  by  equations 
(2.1),  (2.2),  and  (2.3)  has  a  unique  solution.  In  the  following,  we  will  attempt  to  trace  the 
dependence  of  various  constants  such  as:  wave  number  k ,  mesh  size  h,  and  polynomial  order 
of  approximation  p.  C  will  denote  a  generic  positive  constant  independent  of  h,  p,  and  k. 
A  subindex  h,  p,  or  k  will  denote  dependence  upon  h,  p ,  or  k ,  respectively.  For  example, 
Cp  will  denote  a  generic  positive  constant  independent  of  h  and  k,  but  possibly  dependent 
upon  p. 

We  assume  that  our  subspace  decomposition  M  =  J2  Mi  is  such  that  the  discrete  Friedrichs 
inequality  holds,  i.e.: 

||  q i  ||l2<  Ch  ||  Vxq(  ||L2  Vq t  G  Mt  ,  l  >  1  ,  (2.31) 

where  Mt  =  {q;  G  M,  :  (q,,  V(j>i)L 2  =  0  V0;  G  Wi},  and  VIT)  =  {qi  G  M;  :  V  x  qt  =  0}  C  Mt. 
This  inequality  has  been  proved  for  a  variety  of  space  decompositions,  including  spaces 
corresponding  to  local  Dirichlet  problems  for  ftp-meshes  (see,  for  example,  [10,  11,  14]). 

Using  the  discrete  Friedrichs  inequality  we  can  prove  the  following  result. 


Proposition  1 

Let  q^  G  Mi  (/  >1)  be  a  solution  of  problem, 

A(q«,Vi)  =  A(u,  Vvi  G  Mi  ,  (2.32) 

where  u  G  M,  and  A(,  )  is  the  bilinear  form  associated  to  our  variational  formulation.  Then: 


I  qi  I iT£>(curl,f2;)  —  (~fik  ||  U  ||/fD(curl,0;)  ? 
where  supp  q /  C  0/ ,  V  q^  G  M{\,  and  Chk  = 


max{l,  k2} 


(2.33) 

It  follows  that,  if 


min{l  —  C'2h2(k2  +  1  ),k2} 
h2(k2  +  1)  is  small  enough,  then  the  local  problems  have  unique  solutions,  and  therefore,  the 
corresponding  block  Jacobi  smoothers  are  well  defined.  H 
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Proof: 

Using  discrete  Helmholtz  decomposition,  we  have: 
q i  =  q  1,1  +  q/g  ,  (qu  ,  qu)u  =  0  (2.34) 

where  q^i  G  Mi,  and  q^2  =  V0  for  some  (f)  G  IF) . 

Substituting  V;  =  q/  in  (2.32),  and  recalling  the  definition  of  bilinear  form  A(  ,  ),  we 
obtain: 

(Vxqu,  V xqu)i2  -  k2(ql:1,qitl)L2  -  k2( qZ)2,  q/, 2)1,2  = 

(2.35) 

(V xu,  Vxqu)L2  -  k2 (u,  qu)L2  -  k2( u,  q;, 2)l*- 
Setting  v;  =  qi  2  in  (2.32),  we  have: 

k2(qi;2,qia)L^  =  /c2(u,q/i2)L2  .  (2.36) 

Thus: 

(qu,  qi,i)tfD(curi)  -  ( k 2  +  l)(qz,i,  q/,i)i2  =  (V xu,  V xq;>1)L2  -  k2( u,  qM)z,2  .  (2.37) 

Since  q^i  is  discrete  divergence  free  (be.,  q^  G  Mi),  we  can  apply  Friedrich’s  inequality  on 
the  left  hand  side: 


II  qu  HL(curi)  -(^2  + 1)  II  qu  \\b>  [1  -  C2h2(k2  +  1)]  ||  qu  ||^(curl)  .  (2.38) 

Dividing  last  equation  by  ||  q/:)  ||  HD(CUri),  and  applying  eq.  (2.37),  we  obtain: 

h  ^2/202,^111  I,  ^  (Vxu,  Vxqp)p  —  k2(u,  qu)z,2 

[1  C  h  (k  +  1)J  ||  q^i  1 1 Hd (curl) ^  sup  n  n  *  (2.39) 

q  1,16*!  II  ^,1  II Hd(  curl) 

From  eq.  (2.36),  we  derive  for  q/)2  that: 

/2  ii  n  /  (Vxu,  Vxqi  2)L2  —  k2(u,  q/)2)L2 

k  II  q/,2  ||tfD(curl)<  sup  - n - 77 -  • 

q/,2 =vy, yeW;  ||  q/,2  ||Ho(curi) 

Using  (2.39),  (2.40),  and  the  orthogonality  of  Hilbert  spaces  Mi  and  VH7/,  we  conclude 
[1  -  C2h2{k 2  +  l)]2  ||  q/,i  ||^(curl)  +kA  ||  q/,2  ||^(curl)< 


(2.40) 


(Vxu,  Vxq/)i2  -  A;2(u,q/)L2  2  2  2 

(sup  - 77 - 77 - )  <  (||  V  X  U  ||L2  +k  ||  U  ||L2)  , 

q/SM;  II  q/  1 1  Hu  (curl) 


(2.41) 
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and  the  result  follows,  since: 


n 2  t  n  ||2  \i/2  ^  max{l ,k2} 

q/’i  llff°(curl)  +  11  q'’2  H^lcuri))  <  m[n{i^c2h2(k2  +  l),k2} 


I 


(2.42) 


As  a  consequence  of  the  proposition,  the  A(,  )-projection  P/  :  //^(curl)  — >  M/  satisfies: 


max{l,&2} 

P/U  ^D(curl)-  min{l  —  C2h2{k2  +  1),  k2}  11  U  ^D(curI) 


(2.43) 


At  this  point,  we  consider  an  auxiliary  boundary  value  problem  given  again  by  equations 
(2.1),  (2.2),  and  (2.3),  but  with  k 2  =  —1.  And  we  denote  operators  associated  to  our  SPD 
auxiliary  problem  with  the  ~  symbol.  For  example,  P/. 

In  the  following,  we  show  that  the  convergence  properties  of  the  two  grid  solver  for  the 
original  and  auxiliary  problems  are  comparable  up  to  a  perturbation  term.  Such  results  were 
first  proved  for  the  Helmholtz  equation  in  [6].  That  they  can  be  extended  to  the  Maxwell 
case,  notwithstanding  the  nonellipticity,  was  realized  in  [15].  The  following  perturbation 
lemma  is  proved  along  the  lines  of  a  similar  result  in  [16]: 


Lemma  1 

For  all  l  >1,  we  have: 


where  Chk  = 


P l  —  P l  || ipjj (curl) F  C(  1  +  Chk)(k2  +  1  )h  , 
max{l,  k2} 


min{l  —  C2h2(k2  +  1),  k2} 


I 


(2.44) 


Proof: 

Let  u,  q  G  //^(curl).  Since  Pi  is  an  /J/)(curl)-projection,  we  obtain  the  following 
identity: 

(P/U  -  PiU,  q) ipp (curl)  =  (P/U  -  u,  P/q)tfD(curl) 

=  A(Piu  -  u,  P/q)  +  (k2  +  l)(PiU  -  u,  P/q)L2  (2.45) 

=  {k2  +  l)(P/u  -  u,  P;q)L2 

Now,  using  discrete  Helmholtz  decomposition  (and  notation  of  Proposition  1): 

(u  -  P/u,  q/)L2  =  (u  -  P/u,  q/,i)L2  +  (u  -  P/u,  q/,2)n2  .  (2.46) 
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Since  q/j2  =  V^: 


~k2( u  -  P*u,  q;,2)L2  =  A(u  -  P/U,  qij2)  =  0  .  (2.47) 

Applying  Cauchy-Schwarz  inequality,  followed  by  discrete  Friedrich’s  inequality,  we  obtain: 
(u-P,u,qM)L2  <  Ch  ||  u-P,u  ||L2(nj)||  V  x  qM  ||L2(np  .  (2.48) 

Thus: 


(P/U  -  P*u,  q)HD(curi)  <  Ch{k2  +  1)  II  u  P;U  ||L2(ni)||  V  x  P,q  ||L2(n;)  •  (2.49) 

From  (2.43)  we  obtain: 

II  u  —  P/u  IU2(ni)<  (1  +  Chk )  ||  u  ||ffD(curiA)  •  (2.50) 


Finally, 


II  v  x  IU2(ni)<||  q  II hd(  curl,f2/)  •  (2.51) 

And  the  result  follows.  1 

The  following  lemma  quantifies  the  difference  between  the  definite  and  the  indefinite 
coarse  solves  that  appears  in  the  algorithm.  The  proof  is  similar  to  the  proof  of  [15, 
Lemma  4.3]  (also  cf.  [23])  but  we  now  keep  track  of  the  dependence  of  the  constants  on 
polynomial  degree  using  the  recent  results  in  [10,  11]. 


Lemma  2 

If  domain  is  convex,  we  have: 


Po  -  Po  || hd(  curl)  <  Ck 


h 


pL 


v2-*j i  -  a 


h 


kv  l/2-e 


(2.52) 


where  e  >  0.  I 

Proof:  Equations  (2.45)-(2.47)  are  valid  also  for  coarse  grid  subspace  M0.  We  have  for 

all  u,  q  G  M  : 

(PqU  -  P0u,  q) #£> (curl)  =  ( k 2  +  l)(u  -  P0u,  -q0,i)L2  ,  (2.53) 
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where  q0ii  is  the  discrete  divergence  free  part  of  q0  =  Poq. 

We  dehne  e  =  P0u  —  u,  with  e1?  e2  denoting  the  corresponding  terms  of  the  discrete 
Helmholtz  decomposition  of  e.  Then: 

(P0U  —  P oU,  q) Hq (curl)  =  (k2  +  l)[(ei,  q0ii)i,2  +  (e2,  qo,i)i2]  (2-54) 

In  order  to  estimate  both  terms  on  the  right  hand  side,  we  show  first  (following  [15])  that: 

II  ei  ||l2<  Ck  II  e  \\hd (curl)  •  (2.55) 

We  dehne  e^q0,1  G  //^(curl)  by  the  following  conditions: 

V  x  e1  =  V  x  ei  .  ;  (e1  ,  V0)  =  0  V0  G  Hf,  . 

(2.56) 

V  X  q0’1  =  V  X  q0,i  .  ;  (q0'1  ,  V0)  =  0  V0  G  Hf  . 

The  following  result  has  been  proved  in  [3]  for  square  elements  using  the  technique  of 
projection-based  interpolation1 , 

II  el  —  el  |U2<  ^ pl/2-e  II  ^  X  el  IU2  •  (2-57) 

Here  e  >  0  is  an  arbitrary  small  number,  and  C  =  (7(e)  — »  oo,  as  e  — »•  0. 

With  discrete  divergence  free  e,  replaced  with  pointwise  divergence  free  held  e1,  we  turn 
now  to  the  standard  duality  argument  and  consider  solution  w1  G  //^(curl)  to  the  dual 
problem: 

H(p,  w1)  =  (e1,  p)  V  p  G  HD(cur\)  .  (2.58) 

With  the  assumption  on  convexity  of  the  domain,  we  have: 

II  W  | Ih1  (curl)  —  Ck  ||  ®  ||l2  j  (2.59) 

where  ||  w1  ||^i(curl)=||  w1  ||^  +  ||  V  x  w1  ||^. 

Using  the  standard  duality  argument,  and  the  fact  that  w 1  is  divergence  free,  we  obtain, 

II  e1  || l2<  ^(e1,  w1)  =  H(ei,  w1)  =  H(e.  w1)  =  H(e,  w1  —  n^w1) 

(2.60) 

<  Ck  ||  e!  ||/fD(curl)  •  II  W1  -  n^'w1  ||ffD(curl)  , 

1An  analogous  result  holds  for  triangular  elements  under  a  conjecture  of  an  instability  result,  see  [4] 
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where  \[curl  is  the  projection-based  interpolation  operator. 

From  the  approximation  theory  in  [10,  11],  we  obtain: 

||  W1  _  nC“Ww'  || .ffp (curl)  —  C^r~e  II  wl  Hi?1  (curl)  —  Ck^Te  II  el  IU2  ■  (2-61) 

Application  of  triangle  inequality  finishes  proof  of  (2.55). 

Now,  since  q0,1  is  divergence  free,  using  Cauchy-Schwarz  inequality,  and  estimate  (2.57) 
for  function  q0,1,  we  obtain: 


(e2,qo,i )l2  =  (e2,qo.i  -q0,1)i2  <||  e2  lUHl  qo,i  -  q0,1  IU2 


<  c- 


h 


1/2— e 


V 


e  IIl2 II  ^  x  qo,i  IU2<  C- 


h 


1/2— e 


P 


II  Hd{  cur  oil  q  I \hd{  curl) 


(2.62) 


Similarly,  using  Cauchy-Schwarz  inequality,  and  estimate  (2.55),  we  get: 


(el,q0,l)i2  <||  ei  ||i2||  qo.l  |U2<  Ck  _1/2-e  II  e  1 1  (curl)  1 1  q  ||ffD(curl)  • 


p 


(2.63) 


Thus: 


(P0U  -  P 0U,  q)/fD(curl)  <  Ck 


h 


1/2— e 


P 


II  HD(  cur  oil  q  II hd{  curl) 


(2.64) 


In  order  to  finish  the  proof  for  this  lemma,  it  only  remains  to  show  that  ||  e 
1 

Ck  ,  —  ||  u  ||  Hnicurin  This  can  be  done  as  follows: 

A  -  Clh/pW-' 


I  //D(curl)< 


II  e  HL(curl)  =  ll  U  -  P0U  1 1 #£} (curl) 

=  A(u  —  Pou,  u  —  P0u)  +  (1  +  k2)  (u  —  Poll,  u  —  Pou) L2 

=  (u,  u  -  P0u)#D(curl)  +  (1  +  k2)( u  -  P0u.  P0u  -  Pnu)L2  ('2‘65') 

—  II  U  llj?£>(curl)  +Ck  !/2-e  II  U  _  Pou  ||#D(curi) 

I 

Introducing  the  error  reduction  operator  E"  (at  step  n)  associated  to  the  two  grid  algo¬ 
rithm,  he.,  e(n+1)  =  Ene^,  we  conclude  the  following  theorem: 
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THEOREM  1 

If  the  coarse  grid  is  fine  enough,  and 


1 1  E  VI  1 1 Hp (curl) ^  8  ||  U  1 1 Hjj (curl)  with  0  <  8  <  1  ,  (2.66) 

then: 

||  E  u  || //D(curi)  —  8  ||  u  || Hjy(curi)  withO  <  (5  <  1  ,  (2.67) 

where,  8  —  S  +  C max  ||  P/  —  P;  ||ffD(CUri)-  I 

The  remainder  of  this  section  is  devoted  to  proving  (2.66),  which  is  a  sufficient  condition 
to  guarantee  the  main  result  (2.67).  At  this  point,  we  have  already  determined  convergence 
properties  of  the  two  grid  solver  with  respect  to  the  wave  number  k. 

Using  standard  domain  decomposition  techniques  for  SPD  problems,  it  is  well  known  that 
(2.66)  follows  from  the  following  two  conditions  for  the  subspace  splitting  (see,  for  instance, 
[29,  5],  or  [27]): 

•  Condition  I,  necessary  to  estimate  the  maximum  eigenvalue  of  the  preconditioned 
matrix: 

E[(qp  q?)-Hn(curl)  ^  C[^(qj,  q.j)/fD(curl)]  /  [E[(qp  9b)  jjp(curl)]  ^  ;  (2.68) 

i> 1 j> 1  i> 1  j>  1 

where  q,;  G  Mi.  This  condition  is  easily  proved  by  using  a  coloring  algorithm  as 
described,  for  example,  in  [27]. 

•  Condition  II.  necessary  to  estimate  the  minimum  eigenvalue  of  the  preconditioned 
matrix,  by  showing  that  M  =  Jf,  Mi  is  a  stable  subspace  splitting,  he.,  for  all  q  G  M, 
there  exist  q^  G  M;  such  that  q  =  J2i>0  qo  and 

E  II  II  (curl)  —  C  II  Q  1 1 (curl)  ■  (2.69) 

l>  0 

In  order  to  prove  Condition  II  with  a  constant  independent  of  h  (but  possibly  dependent 
upon  p),  we  consider  again  the  discrete  Helmholtz  decomposition.  For  all  q  G  (I  —  Pq )M: 


q  =  'Vw  +  q 


(2.70) 


where  q  G  M  =  {q  G  M  :  (q,  V0) L2  =  0  V0  G  W},  w  G  W . 
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If  domain  II  is  convex,  it  is  proved  in  [2]  using  the  fact  that  q  G  (I  —  P0)M  that: 

I  q  |l2<  Cph  ||  q  ||irD(curi)  •  >.  ||  w  ||r,2<(  Cph  ||  q  ||^2  .  (2.71) 

Assuming  that  we  have  an  L2-stablc  splitting  for  w  (see  [2]),  and  using  discrete  Poincare 
and  inverse  inequalities,  we  obtain: 

II  q  ||l2>  Cph II  w  \\L2>  Y Cph'1  II  Wi  ||L2>  YCP  II  IU2  •  (2-72) 

i>i  i>i 

Similarly,  assuming  that  we  have  an  L2-stable  splitting  for  q,  and  using  discrete  Friedrichs 
and  inverse  inequalities,  we  obtain: 

||  q  |kD(curi)>  Cpil  +  h-1)  ||  q  ||L2>  (2.73) 

Ycp(l  +  h-1)  ||  qi  ||l2>  Y^p  II  ^  II (curl)  •  (2.74) 

i>i  i>i 

Defining  q^  =  q/  +  'Vivi,  we  conclude  that  for  all  q  G  (I  —  Po )M  there  exists  a  decomposition 
q  =  Yj[>i  q i  such  that: 

II  q  II  (curl)  —  CpY  II  ^  1 1  (curl)  ■  (2-75) 

;>i 

Then  Condition  II  holds. 

Observation:  We  have  shown  that  a  sufficient  condition  for  convergence  of  the  two  grid 
solver  is  to  have  an  L2-stable  subspace  splittings  for  both  parts  of  the  discrete  Helmholtz 
decomposition.  In  particular: 

•  Mi  =  Mi  (as  defined  in  Section  2.2),  implicitly  generate  L2-stable  splittings  for  the 
discrete  divergence  free  and  the  gradient  parts. 

•  Mi  =  Mf  for  1  <  l  <  Ne,  Mi  =  VI'F/''  for  Ne  + 1  <  l  <  Ne  +  Nv  (as  dehned  in  Section 
2.2),  generate  L2-stable  splittings  for  the  discrete  divergence  free  (by  using  the  hrst 
Ne  subspaces)  and  the  gradient  parts  (by  using  the  last  Nv  subspaces).  Notice  that 
if  the  last  Nv  subspaces  are  not  included,  then  we  do  not  obtain  a  stable  splitting  for 
gradients,  leading  to  a  diverging  two  grid  algorithm. 

•  The  subspace  splitting  corresponding  to  the  definition  of  our  smoother  Se  generates 
a  L2-stable  splitting  for  the  discrete  divergence  free  part,  while  the  subspace  splitting 
corresponding  to  the  definition  of  our  smoother  Sy  generates  a  L2-stable  splitting  for 
the  gradient  part.  Thus,  we  obtain  a  converging  two  grid  solver. 
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Finally,  notice  that  in  order  to  trace  the  dependence  of  constants  upon  p,  we  cannot  use 
inverse  inequalities.  Thus,  this  part  of  the  proof  of  convergence  becomes  rather  challenging 
and  we  have  not  attempted  it.  Numerical  results  indicate  that  dependence  of  the  two  grid 
solver  contraction  constant  upon  p  is,  at  most,  logarithmic. 

2.5  Stopping  criterion 

Our  ideal  stopping  criterion  translates  into  condition, 

\\eW\\B=\\A-1rW\\B<eTOL.  (2.76) 

Obviously,  this  quantity  is  not  computable,  and  a  stopping  criterion  can  only  be  based  on 
an  approximation  to  it. 

3  Implementation 

In  [25],  we  discussed  several  implementation  aspects,  such  as  assembling,  sparse  storage  pat¬ 
tern,  selection  of  blocks  for  the  block-Jacobi  smoother,  and  construction  of  the  prolongation 
(restriction)  operator  for  elliptic  problems.  The  first  two  implementation  aspects  (assem¬ 
bling  and  sparse  storage  pattern)  are  problem  independent,  while  construction  of  elliptic 
operators  (stiffness  matrix,  block-Jacobi  smoother,  and  prolongation/restriction  operators) 
can  be  naturally  extended  to  electromagnetic  problems  by  using  H( curl)  degrees  of  freedom 
(d.o.f.)  instead  of  H 1  d.o.f..  Thus,  most  implementation  details  discussed  in  [25]  remain 
valid  for  EM  problems  as  well. 

In  this  paper,  we  discuss  the  implementation  of  a  new  embedding  operator  from  gradients 
of  H1  into  il(curl)  arising  for  EM  problems.  This  operator  is  needed  to  construct  the  block- 
Jacobi  smoother  for  gradients. 

3.1  The  transfer  matrix  between  H 1  and  H{ curl) 

We  present  now  shortly  the  main  steps  of  the  algorithm  to  construct  the  transfer  matrix  cor¬ 
responding  to  the  embedding  Vff1  C  iJ(curl).  More  precisely,  if  W  C  H 1  and  M  C  iJ(curl) 
denote  the  FE  spaces  with  the  corresponding  FE  basis  functions  given  by  {ei,  ...,er}  <E  W, 
(gi,  -,gs}  e  M,  we  have: 

r  s 

w  =  J2wiei  an(l  E  =  (3.77) 

*= i  j= i 
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We  seek  a  global  matrix  Tl3  such  that 


Vei  =  '£Tjigj,  (3.78) 

3= 1 

which  implies  the  corresponding  relation  between  the  H1-  and  H( curl)-  degrees  of  freedom 
(d.o.f.), 

r 

VEj  =  Y,Tjm-  (3-79) 

i= 1 

The  following  algorithm  exploits  the  technology  of  constrained  approximation  and  general¬ 
ized  connectivities.  For  an  element  K,  the  global  basis  functions  et  and  g 3  are  related  to  the 
element  shape  functions  (f>k  and  ipi, 

e-i\x  —  y ^Cjk<j>k  ,  g?lr  =  Y.  Dji^i  i  (3.80) 

k  i 

where  Cik  and  Dji  are  the  coefficients  corresponding  to  generalized  connectivities  related  to 
irregular  nodes  and  the  constrained  approximation.  Formulas  (3.80)  imply  the  corresponding 
relations  between  local  and  global  d.o.f. 

The  element  transfer  matrix  Ski  relates  element  shape  functions  (j)k  and  i[) i  according  to 
the  formula 

V0fc  =  E«|.  (3-81) 

i 

For  the  parametric  elements  forming  the  de  Rham  diagram,  the  element  transfer  matrix 
is  independent  of  the  element,  and  it  can  be  precomputed  for  the  master  element  shape 
functions. 

Finally,  due  to  the  hierarchical  construction  of  the  shape  functions,  the  master  element 
transfer  matrix  can  be  precomputed  for  the  maximum  order  of  approximation  with  the  actual 
element  matrix  extracted  from  the  precomputed  one.  We  use  a  simple  collocation  method  to 
precompute  Ski- 

For  an  element  K,  we  have: 

Ve*  =  ]T  Tjigj  =  J2  Tji  J2  Dji'ipj,  ,  (3.82) 

3  3  l 

and,  at  the  same  time, 

Ve?  =  y —  yCikY^Si^i  ,  (3.83) 

k  k  i 
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which  implies, 


£2*  (3.84) 

j  l  k  l 

or 

DtT  =  STC  (3.85) 

In  practice,  it  is  not  necessary  to  invert  matrix  D1 .  This  is  clue  to  the  fact  that  for  each  global 
basis  function  gj,  there  exists  at  least  one  element  K .  for  which  restriction  g  j  |  K  reduces  to 
one  of  the  element  shape  functions,  possibly  premultiplied  with  (—1)  sign  factor.  In  other 
words,  in  the  corresponding  row  in  the  matrix  Djl}  there  is  only  one  non-zero  entry. 

The  formal  algorithm  looks  as  follows. 

•  Initiate  TJt  =  0. 

•  For  each  element  K  in  the  mesh: 

—  For  each  local  H( curl)  d.o.f.: 

*  Exit  the  cycle  if  the  local  d.o.f.  is  connected  to  more  than  one  global  d.o.f. 

*  Determine  the  connected  global  d.o.f.  j  and  coefficient  Dji. 

*  For  each  local  H 1  d.o.f.  k: 

■  For  each  connected  global  H 1  d.o.f.  i  set  Tjj  =  0. 

*  End  of  loop  through  local  H 1  d.o.f. 

*  For  each  local  H 1  d.o.f.  k : 

■  For  each  connected  global  Hl  d.o.f.  i  accumulate  T.J%  =  TJt  +  DJ^C^S^- 

*  End  of  loop  through  local  H 1  d.o.f. 

—  End  of  loop  through  local  i/(curl)  d.o.f. 

•  End  of  loop  through  elements. 


4  Numerical  Results  Concerning  the  Two  Grid  Solver 

This  section  is  devoted  to  an  experimental  study  of  convergence  and  performance  of  our 
two-grid  solver  for  EM  problems.  The  study  will  be  based  on  three  model  EM  problems  that 
we  will  introduce  shortly.  Using  these  examples,  we  will  address  the  following  issues: 
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•  error  estimation  for  the  two  grid  solver, 

•  the  need  for  controlling  gradients  by  using  Hiptmair  or  AFW  decompositions,  and, 

•  the  possibility  of  guiding  the  optimal  /ip-refinements  with  a  partially  converged  solu¬ 
tion. 

4.1  Examples 

We  shall  work  with  three  EM  examples.  For  each  model  problem,  we  describe  the  geometry, 
governing  equations,  material  coefficients,  and  boundary  conditions.  We  also  display  the 
exact  or  approximate  solution,  and  we  briefly  explain  the  relevance  of  each  problem  in  this 
research. 

4.1.1  Diffraction  of  a  plane  wave  on  a  screen 

We  consider  the  problem  of  a  plane  wave  incident  (at  a  45  degree  angle)  to  a  diffracting 
screen. 

Geometry.  Unit  square  ( [0, 1] 2) -  See  Fig.  1. 


Figure  1:  Geometry  and  solution  (module  of  the  second  component  of  the  electric  held)  of 
the  diffraction  of  plane  wave  on  a  screen. 

Governing  equations.  2D  Maxwell’s  equations. 

Material  coefficients.  Free  space. 

Boundary  conditions.  Dirichlet  boundary  conditions  corresponding  to  the  exact  solu¬ 
tion. 
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Exact  solution.  The  exact  solution  can  be  expressed  in  terms  of  Fresnel  integrals  (see, 
for  example,  [9]). 

H(r,  6)  =  — ^ =  exp7r^4:~jkr{F[\/2kr  sin(0/2  —  7t/8)]  +  F[V2kr  sin(6l/2  +  7t/8)]},  (4.86) 

\/7T 


F(u)  =  {exp-7  71-/4 -y/2[C{^2/iru)  -  j S (y/fyiru)]} ,  (4.87) 

C(z)—jS(z)—  [  exp”1/277-7*2  dt  (Fresnel  Integrals).  (4.88) 

Jo 

Solution  is  displayed  in  Fig.  1. 

Observations.  Solution  of  this  2D  wave  propagation  problem  in  free  space  lives  in 
H( curl),  but  not  in  H1. 


4.1.2  Model  waveguide  example 
Geometry.  See  Fig.  2. 


I 


Figure  2:  Geometry  and  FE  solution  (module  of  second  component  of  the  magnetic  field)  of 
the  model  waveguide  problem. 


Governing  equations.  2D  Maxwell’s  equations. 
Material  coefficients.  Free  space. 
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Boundary  conditions.  Cauchy  boundary  condition  (to  model  the  left  excitation  port), 
absorbing  boundary  condition  (to  model  the  right  port)  and  homogeneous  Dirichlet  BC  (to 
model  the  metallic  top  and  bottom  walls).  See  equations  (5.105)  and  (5.106)  for  details  in 
the  Canchy  BC  formulation. 

Exact  solution.  The  exact  solution  is  unknown.  A  FE  solution  is  displayed  in  Fig.  2. 

Observations.  Solution  of  this  2D  wave  propagation  problem  in  free  space  lives  in 
Ff(curl),  but  not  in  H1.  It  involves  four  singular  reentrant  corners. 

4.1.3  A  3D  Electromagnetics  model  problem 
Geometry.  Unit  cube  ( [0,  l]3) -  See  Fig.  3. 


Figure  3:  Geometry  and  exact  solution  of  the  3D  EM  model  problem  with  k  = 
°3*io«J  (sin(357r/180)  cos(257t/180)ux  +  sin(357r/180)  sin(257r/180)uy  +  cos(357t/180))uz. 

Governing  equations.  Maxwell’s  equations. 

Boundary  conditions.  Dirichlet  at  the  three  faces  adjacent  to  the  origin,  and  Cauchy 
(impedance)  at  the  remaining  three  faces. 

Exact  solution.  The  plane  wave  E{ r)  =  ex P"ik'r>  where 

•  r  =  ux£  +  u yy  +  uz^  is  the  position  vector, 

•  k  =  uxkx  +  ny ky  +  u zkz  with  kx ,  ky,  kz  real  indicates  the  wave  amplitude  and  phase, 
and, 

•  p  =  ux  +  uy  +  uz  determines  polarization  of  the  plane  wave. 
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Observations.  We  use  the  example  to  study  performance  of  the  two  grid  solver  and 
illustrate  necessity  (or  not)  of  large  order  of  approximation  p  for  different  values  of  kx,  ky, 
and  kz. 


4.2  Error  estimation 


In  the  following,  we  focus  on  error  estimation  and  selection  of  the  stopping  criterion  discussed 
in  Section  2.5. 

First,  we  consider: 


aHSEr^  || B 

oWSetW 


Then,  we  define  two  error  estimators: 


En(l) 


a^SEr{n)  ||  B 
a(°'>SEr(°')  ||s 


(4.89) 


(4.90) 


En(  2) 


En(  1) 


1  - 


'g(l)+g(0)  \2 


1  - 


'g(»)+g(n~l)\2 


(4.91) 


where  g(n)  =  JZ-^  • 

Figures  4,  and  5  compare  the  accuracy  of  both  error  estimates  (En(l)  and  En( 2))  for 
different  hp-grids  corresponding  to  the  model  waveguide  and  the  diffraction  of  a  plane  wave 
on  a  screen  problems.  More  numerical  results  comparing  both  error  estimators  can  be  found 
in  [25].  These  results  indicate  that  En{2)  is  a  more  accurate  error  estimator  than  En(  1)  in 
most  (but  not  all)  cases. 

4.3  The  need  for  controlling  gradients  of  potentials 

In  the  following,  we  study  numerically  the  need  for  using  a  subspace  decomposition  for  our 
two  grid  solver  that  provides  control  over  gradients,  either  explicitly  (Hiptmair’s  approach) 
or  implicitly  (the  AFW  approach). 

In  Fig.  6,  we  display  convergence  history  for  the  two  grid  solver  algorithm,  with  and 
without  the  explicit  correction  for  gradients  of  potentials.  If  this  correction  is  not  included, 
we  may  loose  control  over  the  kernel  of  the  curl  operator,  leading  to  a  non  converging  (or 
converging  very  slowly)  two  grid  solver  algorithm.  Notice  that  the  problem  is  induced  by 
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NRDOF=3774 


NRDOF=34161 


Figure  4:  Model  waveguide  example.  Exact  (red)  versus  two  estimated  (cyan  -E(  1)-  and 
blue  -E( 2)-)  errors  for  the  two  grid  solver  with  a  3774  (left)  and  34161  (right)  d.o.f.  mesh. 

NRDOF=13946  NRDOF=45830 


Figure  5:  Diffraction  of  a  plane  wave  on  a  screen  problem.  Exact  (red)  versus  two  estimated 
(cyan  -E(  1)-  and  blue  -E( 2)-)  errors  for  the  two  grid  solver  with  a  13946  (left)  and  45830 
(right)  d.o.f.  mesh. 

the  presence  of  the  curl  operator  in  our  variational  formulation,  and  not  by  the  fact  that 
our  problem  may  be  indefinite. 

In  Figures  7,  and  8  we  display  convergence  history  for  the  two  grid  solver  algorithm 
without  the  explicit  correction  for  gradients  for  a  3D  EM  problem  using  different  smoothers 
Si,  S2,  and  S3,  defined  as: 

•  Si  corresponds  to  AFW  decomposition,  that  is,  a  block  (of  the  Jacobi  smoother) 
corresponds  to  the  span  of  all  basis  functions  whose  supports  are  contained  within  the 
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Figure  6:  Diffraction  of  a  plane  wave  on  a  screen  problem  with  a  purely  imaginary  wave 
number  (thus,  the  problem  is  SPD).  Exact  (red)  versus  two  estimated  (cyan  -E(  1)-  and  blue 
-E( 2)-)  errors  for  the  two  grid  solver  with  a  7471  d.o.f.  mesh  (left  figure).  In  the  right  figure, 
correction  for  gradients  were  not  utilized. 

support  of  a  vertex  basis  functions, 

•  blocks  of  S*2  correspond  to  all  d.o.f.  associated  to  a  particular  (modified)  element,  and 

•  blocks  of  S3  correspond  to  all  d.o.f.  associated  to  all  (modified)  elements  adjacent  to 
a  vertex  node. 

Although  a  block  of  S3  is  larger  than  the  corresponding  block  of  Si,  convergence  of  the 
two  grid  solver  is  only  guaranteed  for  smoother  Si-  Indeed,  only  Si  controls  the  kernel  of 
the  curl  without  the  need  for  an  explicit  gradient  correction. 

Unfortunately,  size  of  patches  associated  to  block  Jacobi  smoother  S\  (corresponding  to 
AFW  decomposition)  are  rather  large  (for  p  »  1).  Thus,  the  corresponding  two  grid  solver 
becomes  quite  expensive,  both,  in  terms  of  memory  requirements,  and  CPU  time. 

4.4  Guiding  ^-adaptivity  with  a  partially  converged  solution 

We  come  now  to  the  most  crucial  question  addressed  in  this  article.  For  EM  problems,  how 
much  can  we  relax  our  convergence  tolerance  for  the  two  grid  solver,  without  loosing  the 
exponential  convergence  in  the  overall  hp-adaptive  strategy?  In  the  remaining  of  this  section, 
we  try  to  reach  a  conclusion  via  numerical  experimentation  only. 

We  work  this  time  with  two  examples:  diffraction  of  a  plane  wave  on  a  screen  problem, 
and  the  model  waveguide  example.  For  each  of  these  problems  we  run  the  /ip-adaptive 
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NRDOF=6084 


NRDOF=6084 


Figure  7:  3D  EM  model  problem.  Exact  (red)  versus  two  estimated  (cyan  -E(  1)-  and  blue 
-E( 2)-)  errors  for  the  two  grid  solver  with  a  6084  d.o.f.  mesh  (without  using  an  explicit 
correction  for  gradients),  for  smoothers  Si  (left  figure)  and  S 2  (right  figure). 


Figure  8:  3D  EM  model  problem.  Exact  (red)  versus  two  estimated  (cyan  -E(l)~  and  blue 
-E( 2)-)  errors  for  the  two  grid  solver  with  a  6084  d.o.f.  mesh  (without  using  an  explicit 
correction  for  gradients),  for  smoother  S3. 

strategy  using  up  to  four  different  strategies  to  solve  the  fine  grid  problem: 

1.  a  direct  (frontal)  solver, 

2.  the  two-grid  solver  with  tolerance  set  to  0.01  (as  described  in  section  4.2), 
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3.  the  two-grid  solver  with  tolerance  set  to  0.1,  and 

4.  the  two-grid  solver  with  tolerance  set  to  0.3. 

As  a  measure  of  the  error,  we  use  an  approximation  of  the  error  in  H (curl)-norm.  We 
simply  compute  the  f/(curl)  norm  of  the  difference  between  the  coarse  and  (partially  con¬ 
verged)  fine  mesh  solution.  The  error  is  reported  relative  to  the  norm  of  the  fine  grid  solution, 
in  percent. 

Finally,  for  each  of  the  cases  under  study,  we  report  the  number  of  the  two-grid  iterations 
necessary  to  achieve  the  required  tolerance. 


10" 


512 


Figure  9:  Diffraction  of  a  plane  wave  on  a  screen  problem.  Guiding  hp-refinements  with  a 
partially  converged  solution.  The  left  figure  displays  a  discretization  error  estimate.  The 
right  figure  shows  the  number  of  iterations  needed  by  the  two  grid  solver. 

From  Figures  9  and  10  we  draw  the  following  conclusions. 

•  The  two-grid  solver  with  0.01  error  tolerance  generates  a  sequence  of  hp-grids  that  has 
similar  convergence  rates  to  the  sequence  of  /rp-grids  obtained  by  using  a  direct  solver. 

•  As  we  increase  the  two-grid  solver  error  tolerance  up  to  0.3,  the  convergence  rates  of 
the  corresponding  sequence  of  hp- grids  does  not  decrease  at  all.  However,  the  number 
of  iterations  (as  well  as  the  CPU  time)  decreases  dramatically. 

•  The  number  of  iterations  required  by  our  two-grid  solver  does  not  increase  as  the 
number  of  degrees  of  freedom  increases. 
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Figure  10:  Model  waveguide  example.  Guiding  /ip-refinements  with  a  partially  converged 
solution.  The  left  figure  displays  a  discretization  error  estimate.  The  right  figure  shows  the 
number  of  iterations  needed  by  the  two  grid  solver. 

To  summarize,  it  looks  safe  to  relax  the  error  tolerance  to  0.1  (or  even  larger)  value, 
without  loosing  the  exponential  convergence  rates  of  the  overall  hp  mesh  optimization  pro¬ 
cedure. 


5  Electromagnetic  Applications 

We  conclude  our  work  with  a  number  of  more  realistic  EM  examples  related  to  applications  in 
the  area  of  Petroleum  Engineering,  a  simulation  of  a  waveguide  filter  with  six  inductive  irises, 
and  a  dispersion  error  study  in  3D,  critical  for  radar  simulations.  We  focus  on  advantages  and 
limitations  of  our  numerical  technique  combining  the  fully  automatic  /rp-adaptive  strategy 
with  the  two  grid  solver. 

5.1  Modeling  of  Logging  While  Drilling  (LWD)  EM  measuring 
devices 

In  this  section,  we  consider  two  problems  posed  by  the  oil  company  Baker- Atlas2  in  the  area  of 
LWD  EM  measuring  devices:  an  electrostatic  edge  singularity  problem,  and  an  axisymmetric 
battery  antenna  problem. 

2Baker-Atlas,  a  division  of  Baker-Hughes 
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5.1.1  An  electrostatic  edge  singularity  problem 

A  number  of  LWD  instruments  incorporate  EM  antennas  covered  by  metals  with  plastic 
apertures.  Thus,  edge  singularities  for  the  electric  held  may  occur  on  the  boundary  between 
plastic  and  metal. 

As  a  result,  to  find  the  electric  held  near  edge  singularities  may  become  essential.  In 
addition,  edge  singularities  for  the  electric  (or  magnetic)  held  may  occur  in  the  geologi¬ 
cal  formation.  Strength  of  an  edge  singularity  is  dependent  upon  geometry  and  sources. 
The  fully  automatic  /ip-adaptive  algorithm  does  not  only  detects  singularities,  but  it  also 
distinguishes  between  singularities  of  different  strength. 

Here,  we  focus  on  edge  singularities  arising  in  electrostatic  problems  and  we  present  high 
accuracy  approximations  for  the  electric  held  by  considering  the  following  problem  with  an 
edge  singularity,  for  which  analytical  solution  is  known: 

Geometry.  See  Fig.  11. 


Figure  11:  Geometry  of  the  electrostatic  problem  with  an  edge  singularity. 

Governing  equation.  Laplace  equation  (—A u  =  0). 

Boundary  conditions,  u  =  —In  r,  where  r  =  ( x 2  +  y 2). 

Exact  solution.  The  analytical  solution  is  known,  and  it  was  provided  to  us  by  Baker- 
Atlas3  (see  eq.  (5.98)). 

Observations.  This  2D  problem  has  a  corner  singularity  located  at  (-1,-1),  which 
corresponds  to  an  edge  singularity  in  3D  located  at  (-1,-1  ,z).  We  are  interested  in  ap¬ 
proximating  the  exact  solution  at  distances  from  the  singularity  ten  to  twelve  orders  of 

3We  thank  Lev  Tabarovsky  and  Alex  Bespalov  for  this  example 
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magnitude  smaller  than  the  size  of  the  domain. 

To  introduce  the  physics  of  our  problem,  we  consider  two  perfect  electric  conductor 
(PEC)  infinite  lines  intersecting  at  a  nonzero  point  (xc,yc),  as  shown  in  Fig.  1.  Let  f3  be 
the  angle  between  PEC\  and  PEC2,  and  p  a  Dirac’s  delta  function  distribution  of  charges 
concentrated  at  point  (0,0). 


Figure  12:  Model  problem. 


Boundary  Value  Problem  (BVP).  The  electrostatic  phenomena  is  governed  by  the 
following  system  of  equations, 


Vx£  =  0 
V  •  (eE)  =  p  . 

Thus,  solving  for  scalar  potential  p,  we  obtain: 

div'Vp  =  A  p  =  -  (5.93) 

with  boundary  conditions: 


•  p  —  0  on  PEC  1  and  PFC2. 

•  p  —  0  on  a  boundary  far  enough  from  the  source.  Since  electric  field  E  decays  as 
1/r  (where  r  is  the  distance  from  a  given  point  to  the  source),  for  r  large  enough  (for 
example,  r  =  106),  the  electric  held  intensity  is  negligible,  and  can  be  replaced  by  0. 


Thus,  denoting  by  our  computational  domain  (shown  in  Fig.  11)  and  by  T  its  boundary, 
we  obtain: 


A p  —  -  in  Q 
e 

p  —  0  on  T  . 


(5.94) 
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Secondary  (scattered)  field.  Now,  let  p  be  a  unitary  electric  charge  distribution  con¬ 
centrated  at  the  origin  in  free  space.  Solution  is  given  by  pPrim  =  In  r.  Then: 


A p  -  pPrim  =  0  in  n 
p  —  pPrim  =  0  —  In  r  =  —In  r  on  T 


(5.95) 


Solving  the  problem  for  secondary  potential  psec  =  p  —  pPrim  ,  we  avoid  modeling  of  the 
source,  obtaining: 


A psec  —  0  in  Q 

psec  _  _in  r  on  p 


(5.96) 


Variational  Formulation.  In  order  to  derive  the  corresponding  variational  formulation, 
we  multiply  equation  A psec  =  0  by  a  test  function  v  e  V  —  Hq(Q)  =  {u  G  i/1(h2)  :  u  = 
0  on  T },  and  integrate  (by  parts)  over  domain  hi.  We  obtain: 


Find  u  G  Uq  +  V 
(  VuVv  =  0  Vu  e  V , 


(5.97) 


where  Uq  is  a  lift  corresponding  to  the  non-homogeneous  Dirichlct  boundary  conditions. 


Exact  Solution.  The  exact  solution  of  this  electrostatic  problem  can  be  computed  analyt¬ 
ically  [22],  At  points  located  on  the  surface  of  PEC i,  the  normal  component  of  the  electric 
field  as  a  function  of  (3, 7,  x,  y ,  xc,  and  yc,  is  given  by: 


2t r 


^sin(^) 


En  = 

s  [1-2 lecos(^)+lT] 

where  s  is  the  distance  from  a  given  point  (x,y)  to  the  corner  (xc,yc),  i.e., 
s  =  yj {x  -  xc )2  +  (y-  yc )2 
and  l  is  given  by: 

l  _  \J% 2c  +  y2c 


In  particular,  for  7  =  f ,  we  have: 


En  = 


2vr 


-2vr 


x2c  +  yl\l  0+l  v]  [(Xc  +  Vc)  20 s  0  +  (x2c  +  yl)ws  p . 


(5.98) 


(5.99) 


(5.100) 


(5.101) 
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There  exists  a  singularity  at  point  (xc,  yc)  (edge  singularity  in  three  dimensions)  if  and 
only  if  7 t  <  (3.  Furthermore,  the  larger  is  [3  the  stronger  is  the  singularity.  In  particular,  the 
strength  of  singularity  is  independent  of  the  selected  nonzero  point  ( xc,yc ).  Therefore,  we 
set  xc  —  yc  —  —  1,  obtaining: 


En  =  En(s,  (3) 


-2vr 

r  7r  1  I  7T  1  7T  • 

[2  2/3S1+/3  _|_  22/3  51  /?] 


(5.102) 


For  simplicity,  we  will  restrict  ourselves  only  to  the  case  (3  =  358  degrees,  which  corre¬ 
sponds  to  a  strong  edge  singularity. 


Numerical  Results.  Figures  13,  14,  15, 
final  hp- grid  (generated  automatically  by 
Notice  that  elements  near  the  singularity  ar 
other  elements  in  the  same  grid! 


16,  17,  18,  and  19,  show  different  zooms  on  the 
our  refinement  strategy)  toward  the  singularity, 
up  to  thirteen  orders  of  magnitude  smaller  than 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


Figure  13:  Final  /ip-grid  (Zooms  =  1,10)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 

Finally,  in  Fig.  20,  a  comparison  between  the  exact  and  approximate  solutions  of  the 
normal  component  of  the  electric  field  over  PEC \  at  points  located  near  the  singularity  is 
displayed.  Notice  that  in  order  to  study  behavior  of  the  edge  singularity,  we  are  interested 
in  points  located  at  distances  10~6  —  10-2  from  the  singularity.  As  it  can  be  concluded  from 
Fig.  20,  we  do  fully  resolve  the  problem  in  the  region  of  interest. 
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2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code  2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


Figure  14:  Final  lip- grid  (Zooms  =  102,103)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 


Figure  15:  Final  hp- grid  (Zooms  =  10  ,105)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 

5.1.2  An  axisymmetric  battery  antenna  problem:  the  need  for  goal-oriented 
adaptivity 

We  consider  the  following  axisymmetric  battery  antenna  in  a  homogeneous  medium  with 
nonzero  conductivity  a. 

Geometry.  3D  axisymmetric  problem.  See  Fig.  21. 
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2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


Figure  16:  Final  hp- grid  (Zooms  =  106,107)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  l  (dark 
blue)  to  p  =  8  (pink). 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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Figure  17:  Final  hp- grid  (Zooms  =  108,109)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 


Governing  equations.  Axisymmetric  3D  Maxwell’s  equations. 

Material  coefficients. 

•  Conductivity:  1  Siemens. 

•  Frequency:  1  Mhz. 
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2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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Figure  18:  Final  hp- grid  (Zooms  =  10lo,10n)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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Figure  19:  Final  hp- grid  (Zooms  =  1012,1013)  for  the  electrostatic  edge  singularity  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  =  1  (dark 
blue)  to  p  =  8  (pink). 

•  Relative  permeability  and  permittivity:  1. 

Boundary  conditions.  Homogeneous  Dirichlet  and  Neumann,  see  (5.103). 

Exact  solution.  The  analytical  solution  is  unknown,  although  it  is  known  to  decay 
exponentially  as  we  go  away  from  the  battery  antenna,  since  we  have  a  nonzero  conductivity. 
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58.199630 


Figure  20:  Value  on  the  PEC\  edge  of  normal  component  of  electric  field  at  distances 
1CT6  —  1CT4  (top  figure),  1CT4  —  1CT2  (bottom  left  figure),  and  1CT2  —  10°  (bottom  right 
figure)  from  the  singularity.  The  blue  curve  denotes  exact  solution,  while  FE  approximation 
is  represented  by  a  black  curve. 


Observations.  We  are  interested  in  approximating  the  exact  solution  at  points  far  away 
(0.5m)  from  the  antenna. 

The  original  3D  problem  can  be  reduced  to  a  2D  boundary  value  problem,  which  can  be 
formulated  in  terms  of  a  2D  E-field  obtained  by  solving  the  reduced  wave  equation  with  the 
appropriate  boundary  conditions. 

Since  a  ^  0,  it  is  known  that  the  electric  field  intensity  decays  exponentially  as  we 
move  away  from  the  source  (battery  antenna),  and  we  can  select  a  finite  computational 
domain  with  homogeneous  Dirichlet  boundary  conditions  enforced  at  points  distant  from 
the  antenna.  More  precisely,  we  may  formulate  our  boundary  value  problem  in  domain 
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Figure  21:  A  2D  cross  section  geometry  of  the  axisymmetric  battery  antenna  problem. 


(shown  in  Fig.  21)  as  follows. 


V  x  xhj-  (u2e 

<  n  x  V  x  E  =  —  jcop 
n  x  V  x  E  =  0 
n  x  E  =  0 


jua)E 


0  in  D 
on  Ti 

on  r3 

on  r2  u  r4 . 


(5.103) 


The  corresponding  variational  formulation  is  given  by: 


Find  E  G  F/£)(curl;  fl)  such  that 

[  -(V  x  E)  •  (V  x  F)dx- 
Jn  /i 


(uj2e  —  juja)E  ■  Fdx  = 


u2fi 


i  ri 


FdS  \  for  all  F  G  ^(curl;  D)  . 


(5.104) 


where  Hn( curl;  D)  =  {E  G  H( curl;  O)  :  nxF  |(r2ur4)=  0}  is  the  Hilbert  space  of  admissible 
solutions. 

The  corresponding  stabilized  variational  formulation  can  then  be  derived  using  techniques 
of  Section  2.1. 

Fig.  22  shows  convergence  history  for  a  sequence  of  optimal  hp-grids  produced  by  our 
refinement  strategy,  which  is  delivering  exponential  convergence  rates. 

The  final  /?,p-grid  (with  the  corresponding  zooms  toward  the  battery  antenna)  is  displayed 
in  Figures  23,  24,25,  26,  and  27.  This  final  hp- grid  is  optimal  in  the  sense  that  it  minimizes 
the  relative  energy  norm  error  with  respect  to  the  number  of  unknowns. 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


90J^fcrror 

71.9  2 

573.5 

45/74 

36.47 

29A>? 

23.20 

1830 

14/75 

11.76 

9.38 


SCALES:  nrdofA0.33,  log(error) 


nrdof 

T83  I  579  I  1325  I  2534  I  4317  I  6786  I  10053  114229 


Figure  22:  Convergence  history  of  the  axisymmetric  battery  antenna  in  the  scales  number 
of  unknowns  to  the  power  of  1/3  (algebraic  scale)  vs  logarithm  of  the  relative  energy  norm 
error. 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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Figure  23:  Final  hp- grid  (Zooms  =  1,10)  for  the  axisymmetric  battery  antenna  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  =  1  (dark 
blue)  to  p  =  8  (pink). 


Unfortunately,  the  objective  of  this  problem  is  to  estimate  the  electric  held  at  distances 
far  away  from  the  source.  And  for  this  purpose,  an  energy  norm  based  refinement  strategy 
is  not  suitable,  and  a  goal-oriented  adaptive  strategy  is  needed. 
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2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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Figure  24:  Final  hp- grid  (Zooms  =  102,103)  for  the  axisymmetric  battery  antenna  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  =  1  (dark 
blue)  to  p  =  8  (pink). 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 
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Figure  25:  Final  hp- grid  (Zooms  =  104,105)  for  the  axisymmetric  battery  antenna  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 

5.2  Waveguide  design 

In  this  section,  we  focus  our  attention  on  the  following  six  inductive  irises  waveguide  problem: 
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2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


2Dhp90:  A  Fully  automatic  hp-adaptive  Finite  Element  code 


Figure  26:  Final  hp- grid  (Zooms  =  106,10')  for  the  axisymmetric  battery  antenna  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  =  1  (dark 
blue)  to  p  =  8  (pink). 


Figure  27:  Final  hp- grid  (Zooms  =  108,109)  for  the  axisymmetric  battery  antenna  problem. 
Different  colors  indicate  different  polynomial  orders  of  a  approximation,  from  p  —  1  (dark 
blue)  to  p  =  8  (pink). 

Geometry:  A  six  inductive  irises  filter4  of  dimensions  «  20  x  2  x  1  cm.  For  details,  see 
Fig.  28. 

Governing  equations.  Maxwell’s  equations. 

Material  coefficients. 

4We  thank  Dr.  Luis  Garcia-Castillo  and  Mr.  Sergio  Llorente  for  this  problem 
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Figure  28:  Geometry  and  FE  solution  (^|  Hx  |2  +  |  Hy  |2)  at  8.82  Ghz  of  the  waveguide 
problem  with  six  inductive  irises. 

•  Free  space 

•  Operating  Frequency  fa  8.8  —  9.6  Ghz 

•  Cutoff  frequency  fa  6.56  Ghz 

Boundary  conditions.  Dirichlet,  Neumann  and  Cauchy,  see  (5.105). 

Exact  solution.  The  exact  solution  is  unknown.  A  FE  solution  is  displayed  in  Fig.  28. 

Observations. 

•  H-plane  six  inductive  irises  filter. 

•  Solution  of  this  wave  propagation  problem  lives  in  //(curl),  but  not  in  H1. 

•  Solution  involves  resolution  of  twenty  four  singular  reentrant  corners. 

•  Dominant  mode  (source):  TE10— mode. 

First,  we  formulate  the  boundary  value  problem.  Next,  we  present  the  variational  for¬ 
mulation  and  discuss  a  way  to  compute  the  scattering  parameters,  the  primary  quantity 
of  interest  for  the  waveguide  design.  Third,  we  display  results  obtained  with  our  numeri¬ 
cal  method.  Finally,  we  use  the  problem  to  illustrate  some  of  the  limitations  of  the  fully 
automatic  hp-adaptive  strategy,  and  the  two  grid  solver. 
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Formulation.  We  excite  the  TEj0-mode  at  the  left  hand  side  port  of  the  waveguide  struc¬ 
ture  (see  Fig.  28),  and  we  consider  an  H-plane  discontinuity,  he.,  the  magnetic  held  is 
invariant  with  respect  to  the  z-direction.  Thus,  our  original  3D  problem  can  be  reduced  to  a 
2D  boundary  value  problem,  which  we  can  formulate  in  terms  of  the  2D  //-field  as  follows. 


V  x  Q  V  xHj-  to2  pH  =  0  in  Q 

n  x  -V  x  H  =  x  n  x  (2 Hinc  -  H)  on  T1 

1  (5-105) 

n  x  -V  x  H  = - — ^ -n  x  n  x  H  on  T2 

|  Pio 

n  x  -V  x  H  =  0  on  r3  , 


v  e 

where  Tx,  T2,  and  T3  are  the  parts  of  the  boundary  corresponding  to  the  excitation  port  (left 
port),  non-excitation  port  (right  port),  and  the  perfect  electric  conductor  respectively  ;  (3 io 
refers  to  the  propagation  constant  of  the  TE10  mode  ;  and  Hinc  is  the  incident  magnetic 
held  at  the  excitation  port. 

The  corresponding  variational  formulation  is  given  by: 


(  Find  H  G  Ho{ curl;  fl)  such  that 


/  -(V  x  H)  ■  (V  x  F)dx  -  /  u2/iH  •  Fdx  + 
In  e  Jn 


ju2H 


/5io  ^riur2 


(n  x  H)  ■  (n  x  F)dS  = 


(5.106) 


(n  x  Hinc)  ■  (n  x  F)dS  for  all  F  G  HD{ curl;  O)  . 


Pio  JTi 

In  the  above  Hd{ curl;  O)  is  the  space  of  admissible  solutions, 


Hd( curl;  Q)  :=  {H  G  L2(Q)  :  V  x  H  G  L2(Q)}  . 


(5.107) 


The  stabilized  variational  formulation  can  be  derived  using  techniques  of  Section  2.1. 


Scattering  parameters.  The  objective  of  the  waveguide  problem  is  to  compute  the  so 
called  scattering  parameters.  Since  the  only  propagating  mode  is  TEW,  we  have  two  power 
waves5  present  at  each  port  Tj  (i  =  1,2):  one  going  inward  (a*),  and  other  going  outward 

5 A  power  wave  can  be  identified  with  a  complex  number  such  that  its  magnitude  squared  represents  the 
power  carried  by  the  wave,  and  its  argument  is  the  phase  of  the  wave 
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(6j)  the  structure.  The  relation  between  the  power  waves  is  linear,  and  may  be  written  in 
matrix  form  as 


h 

b2 


S12  V  (  \ 

S22  )  \  a>2  ) 


(5.108) 


where  Stl  are  the  scattering  parameters,  or  simply  S  parameters.  Thus,  the  S  parameters 
relate  the  incident  and  reflected  power  waves.  Note  that  Su,  S22  are  reflection  coefficients 
and  S'12,5'21  are  transmission  coefficients. 

The  absorbing  boundary  condition  at  T2  implies  a2  =  0,  and  (5.108)  reduces  to: 

S11  =  —  ;  <S'i2  —  —  •  (5.109) 

CL\  CL\ 

Also,  we  have  the  reciprocity  condition,  given  by  S12  =  £21  (see  [17]).  And  since  no  losses 
occur  within  the  waveguide  structure,  symmetry  property  |  An  |J  +  |  S22  |2=  1  holds  (see 
[18]). 


Numerical  results.  The  problem  may  be  solved  by  using  semi-analytical  techniques  (for 
example,  mode  matching  techniques  [28]).  Nevertheless,  it  would  be  desirable  to  solve  it  using 
purely  numerical  techniques,  since  a  numerical  method  allows  for  simulation  of  more  complex 
geometries  and/or  artifacts  possibly  needed  for  the  construction  of  an  actual  waveguide. 

While  attempting  to  solve  this  problem  using  the  fully  automatic  /ip-adaptive  strategy 
coupled  with  the  two  grid  solver,  we  encountered  the  following  limitations. 

1.  We  cannot  guarantee  convergence  of  the  two  grid  solver  if  the  coarse  grid  is 
not  fine  enough,  as  predicted  by  the  theory.  In  this  case,  we  need  a  minimum  of  seven 
elements  per  wavelength  A  in  the  x-direction,  and  thirteen  elements  per  wavelength  A 
in  the  y-direction.  Furthermore,  as  indicated  in  Table  1,  convergence  (or  not)  of  the 
two  grid  solver  is  (almost)  insensitive  to  /^enrichment. 


Does  the  two  grid  solver  converge? 

p  =  1 

p  =  2 

p  =  3 

p  =  4 

Number  of  elements  per  A  =  7, 13 

YES 

YES 

YES 

YES 

Number  of  elements  per  A  =  7, 11 

NO 

NO 

NO 

YES 

Number  of  elements  per  A  =  6, 13 

NO 

NO 

NO 

NO 

Table  1:  Convergence  (or  not)  of  the  two  grid  solver  for  different  quasi-uniform  initial  grids 


2.  We  cannot  guarantee  the  optimality  of  the  fully  automatic  hp-adaptive  strat¬ 
egy  if  the  dispersion  error  is  large.  Since  solution  of  the  problem  on  the  fine  grid 
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is  used  to  guide  optimal  /tp- refinements,  we  need  to  control  the  dispersion  error  on  the 
fine  grid.  Thus,  h  needs  to  be  sufficiently  small  or  p  sufficiently  large.  In  Fig.  29,  we 
compare  the  convergence  history  obtained  by  using  the  fully  automatic  /ip-adaptive 
strategy  starting  with  different  initial  grids.  For  third  order  elements,  the  dispersion 
error  is  under  control  (see  estimates  of  [1,  20,  21]),  and  the  fully  automatic  hp-adaptive 
strategy  converges  exponentially.  We  also  observe  that,  in  the  asymptotic  regime,  all 
curves  present  similar  rates  of  convergence. 
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Figure  29:  Convergence  history  using  the  fully  automatic  /rp-adaptive  strategy  for  different 
initial  grids.  Different  colors  correspond  to  different  initial  orders  of  approximation.  27  is  the 
minimum  number  of  elements  needed  to  reproduce  the  geometry,  while  1620  is  the  minimum 
number  of  elements  needed  to  reproduce  the  geometry  and  to  guarantee  convergence  of  the 
two  grid  solver. 


We  solved  the  six  irises  waveguide  problem  delivering  a  0.2%  error  (in  the  relative  energy 
norm).  Fig.  30  displays  the  magnitude  of  the  Su  scattering  parameter  (on  the  decibel  scale) 
with  respect  to  the  frequency.  This  quantity  is  usually  referred  to  as  the  return  loss  of  the 
waveguide  structure.  For  frequency  interval  8.8  —  9.6  Ghz,  the  return  loss  is  below  — 20dB, 
which  indicates  that  almost  all  energy  passes  through  the  structure,  and  thus,  the  waveguide 
acts  as  a  filter. 

Finally,  Figures  31,  32,  33,  and  34  display  solution  at  different  frequencies.  For  frequencies 
8.72  Ghz,  and  9.71  Ghz,  the  return  loss  of  the  waveguide  structure  is  large,  and  for  frequencies 
8.82  Ghz,  and  9.58  Ghz  the  return  loss  is  below  — 20dB. 
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Figure  30:  Return  loss  of  the  waveguide  structure. 


Figure  31:  |  Hx  |  (upper  figure),  |  Hy  |  (center  figure),  and  ^ \  Hx  |2  +  |  Hy  |2  (lower  figure) 
at  8.72  Ghz  for  the  six  irises  waveguide  problem. 


5.3  Analysis  of  a  3D  EM  model  problem. 

In  this  section,  we  study  numerically  a  3D  EM  model  problem  introduced  in  Section  4.1.3. 
Given  a  unit  cube  geometry,  the  objective  is  to  determine  number  of  elements  N  and  corre¬ 
sponding  order  of  approximation  p,  needed  to  solve  Maxwell’s  equations  numerically,  for  an 
arbitrary  plane  wave  solution  at  different  frequencies  (thus,  wavelengths). 

A  second  motivation  for  this  numerical  experiment  comes  from  studying  the  fully  auto- 
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Figure  32:  |  Hx  |  (upper  figure),  j  Hy  |  (center  figure),  and  \J \  Hx  |2  +  j  Hy  |2  (lower  figure) 
at  8.82  Ghz  for  the  six  irises  waveguide  problem. 


Figure  33:  |  Hx  |  (upper  figure),  |  Hy  |  (center  figure),  and  ^ \  Hx  |2  +  |  Hy  |2  (lower  figure) 
at  9.58  Ghz  for  the  six  irises  waveguide  problem. 


rnatic  hp-adaptive  strategy,  which  may  produce  misleading  results  if  dispersion  error  on  the 
initial  fine  grid  is  too  large.  In  this  section,  we  present  combinations  of  uniform  /rp-grids, 
that  lead  to  solution  of  our  model  problem  (at  different  frequencies)  with  a  relative  energy 
norm  error  below  5%.  These  hp-meshes  may  be  utilized  as  a  guide  to  construct  the  initial 
fine  grid  for  the  /?.p-adaptive  strategy. 

For  these  purposes,  we  select  an  incoming  plane  wave  with 

k  =  /c(sin(a:)  cos(/3)ux  +  sin(a:)  sin(/d)uy  +  cos(/3)uz)  ,  (5.110) 
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Figure  34:  |  Hx  |  (upper  figure),  |  Hy  |  (center  figure),  and  ^ \  Hx  |2  +  |  Hy  |2  (lower  figure) 
at  9.71  Ghz  for  the  six  irises  waveguide  problem. 

where  k  is  the  wave  number,  a  =  357t/180,  and  (3  =  257t/180.  We  consider  a  ID  cross 
section  given  by  the  main  diagonal  of  the  unit  cube,  starting  at  the  origin  and  ending  at 
point  (1,1,1).  Figures  35,  and  36  show  a  comparison  between  exact  and  FE  solution  for 
the  real  part  of  the  first  component  of  the  electric  field,  using  different  hp- meshes.  For 
an  hp- grid  delivering  4.4%  relative  error  in  the  energy  norm,  differences  in  both  phase  and 
amplitude  between  exact  and  FE  solutions  cannot  be  appreciated.  As  the  error  increases, 
these  differences  become  larger. 

Since  solution  of  the  3D  EM  problem  is  smooth,  large  elements  with  large  polynomial 
order  of  approximation  p  are  preferred  over  small  elements  with  small  p.  In  addition,  dis¬ 
persion  error  decreases  faster  by  increasing  p  rather  than  by  decreasing  element  size  h  (see 
[1] , [20] ,  and  [21]).  Table  2  illustrates  these  assertions.  Using  p  =  5,  a  grid  with  51000  d.o.f 
delivers  smaller  error  than  a  grid  with  20  million  unknowns  and  lowest  order  elements  for 
the  model  problem  with  8  wavelengths  on  the  main  diagonal. 

Table  2  has  been  generated  by  using  a  direct  (frontal)  solver.  Notice  that  the  two  grid 
solver  requires  an  elevated  number  of  elements  per  wavelength  on  the  coarse  grid  to  guarantee 
convergence,  as  mentioned  in  Section  5.2.  As  a  consequence,  the  two  grid  solver  could  not 
be  utilized  to  produce  the  table. 
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Nr.  of  A  vs  p 

p—  1 

p— 2 

p= 3 

p— 4 

p= 5 

A=  1 

ERROR 

ELEM/A 

D.O.F. 

5.0  % 

20 

40K 

4.2  % 

3 

946 

1.2  % 

2 

1033 

1.8  % 

1 

308 

0.3  % 

1 

548 

A=  2 

ERROR 

ELEM/A 

D.O.F. 

(>300K) 

4.2  % 

3 

6427 

2.9  % 

1.5 

2764 

1.9  % 

1 

2226 

0.3  % 

1 

4109 

A=  3 

ERROR 

ELEM/A 

D.O.F. 

(>1200K) 

3.5  % 
3.33 

31K 

4.1  % 
1.33 

7115 

1.8  % 

1 

6148 

2.1  % 
0.66 

4109 

A=  4 

ERROR 

ELEM/A 

D.O.F. 

(>2300K) 

(>82K) 

5.0  % 
1.25 

12K 

1.9  % 

1 

14K 

1.2  % 
0.75 
12K 

A=  5 

ERROR 

ELEM/A 

D.O.F. 

(>135K) 

3.6  % 

1.4 

31K 

4.4  % 
0.8 

27K 

3.4  % 

0.6 

12K 

A=  6 

ERROR 

ELEM/A 

D.O.F. 

(>240K) 

4.2  % 
1.33 

46K 

3.8  % 
0.83 

27K 

2.1  % 
0.66 

27K 

A=  7 

ERROR 

ELEM/A 

D.O.F. 

(>410K) 

(>98K) 

3.4  % 
0.86 
45K 

4.3  % 
0.57 
27K 

A=  8 

ERROR 

ELEM/A 

D.O.F. 

(>20M) 

(>650K) 

(>167K) 

(>71K) 

2.8  % 
0.625 
51K 

A=  50 

ERROR 

ELEM/A 

D.O.F. 

(>5000M) 

(>122M) 

(>25M) 

(>14M) 

(>9.5M) 

Table  2:  3D  electromagnetics  model  problem.  For  frequencies  from  1  to  50  wavelengths 
A,  and  uniform  hp-grids  (1  <  p  <  5),  we  display  relative  error  in  the  energy  norm  (in 
percentage),  number  of  elements  per  wavelength  A,  and  actual  (or  estimated)  number  of 
d.o.f  required  to  obtain  an  error  below  5%. 
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RELATIVE  ERROR  IN  THE  ENERGY  NORM=  4.4% 


DISTANCE  FROM  THE  ORIGIN 


RELATIVE  ERROR  IN  THE  ENERGY  NORM=  13.4% 


DISTANCE  FROM  THE  ORIGIN 


Figure  35:  3D  EM  model  problem.  ID  cross  section  over  the  main  diagonal  of  the  unit  cube 
(from  (0,0,0)  to  (1, 1, 1)).  Comparison  between  the  exact  and  the  FE  solution  component 
Re(Ei ),  for  different  hp- meshes  delivering  4.4%  and  13.4%  error  (in  the  relative  energy  norm) 
respectively. 


RELATIVE  ERROR  IN  THE  ENERGY  NORM=  32.8% 


Figure  36:  3D  EM  model  problem.  ID  cross  section  over  the  main  diagonal  of  the  unit  cube 
(from  (0,0,0)  to  (1, 1, 1)).  Comparison  between  the  exact  and  the  FE  solution  component 
Re(Ei ),  for  an  hp- mesh  delivering  32.8%  error  (in  the  relative  energy  norm). 

6  Conclusions  and  Future  Work 

In  this  paper,  we  have  studied  a  two  grid  solver  for  solving  linear  systems  resulting  from 
hp  FE  discretizations  of  Maxwell’s  equations.  The  meshes  come  in  pairs,  consisting  of  a 
coarse  mesh  and  the  corresponding  hne  mesh  obtained  via  the  global  hp  refinement  of  the 
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coarse  mesh.  The  coarse  meshes  are  generated  by  a  special  hp-adaptive  algorithm,  based  on 
minimizing  the  projection  based  interpolation  error  of  the  fine  mesh  solution  with  respect  to 
the  next  optimally  refined  coarse  mesh.  The  solver  combines  block  Jacobi  smoothing  with 
an  optimal  relaxation,  with  the  coarse  grid  solve. 

Instead  of  using  the  two  grid  iteration  for  producing  a  preconditioner  for  Conjugate 
Gradient  (CG)  only,  we  chose  to  accelerate  each  smoothing  operation  individually  with  an 
approximation  of  the  Steepest  Descent  (SD)  method,  which  we  interpret  as  determining  the 
optimal  relaxation  parameter. 

Within  the  described  framework,  we  have  studied  several  critical  questions  including  the 
convergence  theory,  implementation  issues,  the  need  of  controlling  gradients,  error  estimation 
for  the  two  grid  solver  and,  first  of  all,  the  possibility  of  guiding  the  hp  strategy  for  EM 
problems  with  only  partially  converged  solution.  As  a  result  of  it,  we  verified  that  a  partially 
converged  solution,  with  a  rather  large  (relative)  error  tolerance  of  0.1,  is  sufficient  to  guide 
the  hp  strategy.  The  corresponding  number  of  two  grid  iterations  stays  then  very  minimal 
at  a  level  below  10  iterations  per  mesh. 

Then,  we  have  applied  our  numerical  technique  (the  fully  automatic  /ip-adaptive  strategy 
coupled  with  the  two  grid  solver)  to  a  number  of  practical  EM  problems.  While  most 
applications  were  solved  with  extreme  accuracy,  we  also  faced  a  number  of  limitations: 

•  The  two  grid  solver  may  not  converge  for  indefinite  problems  if  the  coarse  grid  is  too 
coarse.  Furthermore  (as  shown  in  Lemma  1),  this  condition  over  the  element  size  does 
not  depend  upon  the  polynomial  order  of  approximation  p.  A  multigrid  solver  for  which 
the  constant  in  Lemma  1  decreases  as  p  increases  is  badly  needed  for  wave  propagation 
problems,  when  hp-Finite  Elements  are  used. 

•  An  adaptive  strategy  based  on  minimization  of  the  energy  norm  may  be  inadequate 
for  a  number  of  EM  applications  as,  for  example,  the  axisymmetric  battery  antenna 
problem.  Thus,  an  hp  goal-oriented  adaptive  algorithm  is  needed. 
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